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DERIVATION  AND  APPLICATION  OF  DUAL-SURFACE 

INTEGRAL  EQUATIONS  FOR  THREE- 
DIMENSIONAL,  MULTI-WAVELENGTH  PERFECT 

CONDUCTORS 


1.  INTRODUCTION 

The  numerical  solution  of  dual-surface  integral  equations'  applied  to  three- 
dimensional  (3-D),  multi-wavelength,  perfectly  conducting  bodies  can  be  obtained  with 
readily  available  computers  in  a  central  processing  unit  (CPU)  time  proportional  to 
approximately  (sA)4  In  (sA)  using  the  conjugate  gradient  method2  and  direct  access 
memory  files3  (s  is  the  dimension  of  the  body  and  \  the  wavelength).  Specifically,  for 
a  given  incident  fieid  and  aspect  angle,  the  induced  cuirent  and  far  field  over  4ji 
steradians  of  a  perfectly  conducting  cube  5  wavelengths  on  a  side  is  computed4  in  about 


(Received  for  Publication  5  Oct  1989) 

1  Tobin,  A.R.,  Yaghjian,  A.D.,  and  Bell,  M.M.  (1987)  Surface  integral  equations  for 
multi-wavelength,  arbitrarily  shaped,  perfectly  conducting  bodies.  Digest  of  the 
National  Radio  Science  Meeting  (URSl),  Boulder,  CO,  p.  9. 

1  Sarkar,  T.K.  and  Arvas,  E.  (1985)  On  a  class  of  finite  step  iterative  methods 
(conjugate  directions)  for  the  solution  of  an  operator  equation  arising  in 
electromagnetics,  IEEE  Trans  Antennas  and  Propagat.  A  P-33:  1058-1066. 

Woodworth,  M  B.  (1988)  Large  Matrix  Solution  Techniques  Applied  to  an  Electro¬ 
magnetic  Scattering  Problem.  RADC-TR-88-268.  ADA206917 

4  Cote,  M.G.,  Woodworth,  M  B.,  and  Yaghjian,  A.D.  (1988)  Scattering  from  the  perfectly 
conducting  cube,  IEEE  Trans.  Antennas  Propagat.  A  P-36:  1321-1329. 
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1.5  hours  of  CPU  time  [utilizing  two-fold  symmetry  of  the  cube:  see  Eq.  (34)]  on  a  VAX 
8650  computer  with  a  "32-bit  Linpack  benchmark  performance  rating"  of  1.3  megaflops.5 
The  same  computer  would  take  about  11  hours  of  CPU  time  for  this  5\  cube  if  the  matrix 
of  the  dual-surface  integral  equation  were  solved  using,  instead  of  the  conjugate 
gradient  method,  Gaussian  elimination,  which  requires  a  computer  time  proportional  to 
approximately  (s/A,)6.  [Gaussian  elimination  CPU  time  with  direct  access  files  can  be 
estimated  from  Eq.  (30)  with  N  =  7 5 (s/A.)  ■]  If  the  two-fold  symmetry  of  the  cube  were 
not  used  to  reduce  the  number  of  unknowns  by  a  factor  of  4  these  CPU  times  using  the 
conjugate  gradient  and  Gaussian  elimination  methods  would  increase  by  factors  of 
approximately  22  and  64,  respectively,  that  is,  from  1.5  and  11  hours  to  about  33  and 
700  hours  (30  days)  of  CPU  time  for  a  5A.  scatterer.  This  latter  CPU  time  of  30  days 
confirms  that  scattering  or  radiation  from  arbitrarily  shaped  5A.,  3-D  bodies  cannot  be 
determined  in  a  reasonable  amount  of  computer  time,  using  conventional  Gaussian 
elimination,  by  a  computer  with  a  Linpack  performance  rating  on  the  order  of  one 
megaflop.  It  becomes  necessary  to  use  faster  matrix  solution  schemes,  such  as  the 
conjugate  gradient  iterative  method,  when  the  integral  equations  are  applied  to 
progressively  larger  bodies,  regardless  of  the  speed  of  the  computer. 

Herein,  we  derive  the  dual-surface  electric  and  magnetic-field  integral  equations  for 
3-D  perfectly  electrically  conducting  bodies,  prove  that  they  produce  a  unique  solution 
at  all  real  frequencies,  and  demonstrate  their  applicability  to  multi-wavelength  bodies 
by  solving  the  dual-surface  magnetic-field  integral  equation  for  a  rectangular  scatterer 
using  the  method  of  conjugate  gradients. 

Magnetic-field  surface  integral  equations  for  perfect  conductors  appeared  in  the 
literature  as  early  as  193 16,  and  both  electric  and  magnetic-field  surface  integral 
equations  were  derived  in  Maue’s  definitive  1949  Zeitschrift  Fur  Physik  paper.7 
However,  only  in  the  last  ten  years  or  so  have  digital  computers  become  fast  enough  to 
solve  these  surface  integral  equations  for  arbitrarily  shaped,  3-D,  multi-wavelength 
bodies. 


5  Dongarra,  J.,  Martin,  J.L.,  and  Worlton,  J.  (1987)  Computer  benchmarking:  paths  and 
pitfalls,  IEEE  Spectrum,  24:  38-43. 

6  Murray,  F.H.  (1931)  Conductors  in  an  electromagnetic  field,  Am.  J.  Math.,  53: 

275-288. 

7  Maue,  A.W.  (1949)  On  the  fonnulation  of  a  general  scattering  problem  by  means  of  an 
integral  equation,  Z.  Phys.,  126(7/9):  601-618. 
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Unfortunately,  as  Murray  and  Maue  noted,  the  original  electric  and  magnetic-field 
integral  equations  (EFIE  and  MFIE)  fail  to  produce  a  unique  exterior  solution  at 
frequencies  equal  to  the  resonant  frequencies  of  the  corresponding  interior  cavity. 

Since  the  density  of  cavity  resonant  frequencies  increases  rapidly  beyond  the  fust 
resonant  frequency,  which  occurs  when  the  dimension  of  a  full-bodied  3-D  scatterer 
equals  about  one  wavelength,  the  numerical  solution  of  3-D,  multi-wavelength  bodies  is 
severely  hampered  by  these  spurious  resonances. 

In  Reference  it  was  proven  that  the  original  integral  equations  allow  spurious 
solutions  at  the  cavity  resonant  frequencies  because  at  (and  only  at)  these  frequencies 
the  MFIE  does  not  restrict  the  tangential  electric  field  to  zero  on  the  surface  of  the 

scatterer  and  the  EFIE  does  not  restrict  the  tangential  magnetic  field  to  K  x  n  on  the 

surface  of  the  scatterer.  (K  is  the  surface  current  and  n  the  outward  unit  normal  to  the 
scatterer.  Interestingly,  the  MFIE  result  was  also  proven  in  the  early  paper  by 
Murray. *)  Among  the  alternatives  that  have  been  proposed  for  eliminating  the  spurious 
solutions  from  the  original  integral  equations,  the  combined-field9’  10  ’  11  or 
combined-source 13  ’  14  integral  equation,  and  the  augmented  electric  or  magnetic-field 


Yaghjian,  A.D.  (1981)  Augmented  electric  and  magnetic-field  integral  equations.  Radio 
Science,  16:  987-1001. 

Mitzner,  K.M.  (1968)  Numerical  solution  of  the  exterior  scattering  problem  at 
eigen-frequencies  of  the  interior  problem.  Digest  of  Fall  URSi  Meeting,  Boston,  MA, 
p.  75. 

10  Poggio,  A.J.,  and  Miller,  E.K.  (1973)  Integral  equation  solutions  of  three-dimensional 
scattering  problems,  in  Computer  Techniques  for  Electromagnetics,  edited  by  R.  Mittra, 
159-264,  Pergamon,  New  York. 

"  Mautz,  J.R.  and  Harrington,  R.F.  (1978)  H-field,  E-field,  and  combined-field 
solutions  for  conducting  bodies  of  revolution.  Arch.  Elektron.  Ubertragungstech. 

< Electron .  Commute),  32(4):  157-164. 

Panic,  I.O.  (1965)  On  the  solvability  of  exterior  boundary-value  problems  for  the  wave 
equation  and  for  a  system  of  Maxwell’s  equations,  Uspehi  Mat.  Nauk,  20(1): 

22!  226. 

M  Brakhage,  H.  and  Werner,  P.  (1965)  Uber  das  Dirichletsche  Aussenraumproblem  fur  die 
Heimholtzsche  Schwingungsgleichung,  Arch.  Math.,  16:  325-329. 

14  Mautz,  J.R.  and  Harrington,  R.F.  (1979)  A  combined-source  solution  for  radiation  and 
scattering  from  a  perfectly  conducting  body,  IEEE  Trans.  Antennas  and  Propagation, 

A P-27:  445-454. 
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integral  equation  appear  the  more  generally  applicable  and  effective  in  numerical 
practice.  However,  for  arbitrarily  shaped,  3-D,  multi-wavelength  bodies,  the  combined 
and  augmented  integral  equations  also  have  their  drawbacks.  The  combined-field  and 
combined-source  equations  involve  the  operators  of  both  the  magnetic-field  equation  and 
the  electric-field  equation,  which  takes  considerably  more  programming  ingenuity  and 
computer  time  than  the  original  MFIE  to  achieve  the  same  accuracy  of  solution.  The 
augmented  MFIE  involves  only  the  magnetic-field  operator,  but  the  augmented  integral 
equations  require  a  special  procedure  to  eliminate  all  the  spurious  solutions  from 
bodies  of  revolution. 

Thus,  we  begin  the  integral  equation  solution  for  arbitrarily  shaped,  3-D, 
multi-wavelength  perfect  conductors  with  the  derivation  of  dual-surface  electric  and 
magnetic-field  integral  equations  that  differ  only  slightly  from  the  original  electric 
and  magnetic-field  integral  equations,  yet  eliminate  all  spurious  solutions.  The 
dual-surface  magnetic-field  integral  equation  was  given  in  Reference  1,  but  the 
derivation  and  proof  of  uniqueness  of  the  dual-surface  magnetic  and  electric-field 
integral  equations  have  not  appeared  previously.  Recently,  Toyoda  et  al.15  presented  an 
"extended  integral  equation  formulation"  for  2-D  scatterers  that  used  additional 
surfaces  near  the  surface  of  the  scatterer.  Their  formulation  for  perfect  conductors 
applies  an  extended  integral  equation  to  an  interior  surface  and  requires  the  interior 
surface  to  move  with  frequency  to  maintain  uniqueness  of  solution.  The  derivation  of 
the  dual-surface  magnetic  and  electric-field  integral  equations,  (5a)  and  (5b),  in  the 
following  section  requires  the  introduction  of  a  second  surface  interior  and  parallel  to 
the  surface  of  the  scatterer,  but  the  resulting  integral  equations  have  a  unique 
solution  at  all  real  frequencies  and  are  applied  to  the  single  surface  of  the  scatterer. 


2.  DERIVATION  OF  DUAL-SURFACE  INTEGRAL  EQUATIONS 

A  time  harmonic  [exp  (-itot),  tn  real  and  >  0]  electromagnetic  field  (E^,  H^c) 
incident  in  free  space  upon  the  surface  S  of  a  perfectly  electrically  conducting 

15  Toyoda,  I.,  Matsuhara,  M.,  and  Kumagai,  N.  (1988)  Extended  integral  equation 
formulation  for  scattering  problems  from  a  cylindrical  scatterer,  IEEE  Trans. 
Antennas  and  Propagat.  36:  1580-1586. 
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scatterer  excites  a  surface  current  K.  (Let  S  be  coincident  with  the  surface  current 

K.)  Since  the  total  field  inside  the  scatterer  is  zero,  the  scattered  fields  equal  the 
negative  of  the  incident  fields  inside  S,  and  one  can  write  the  "interior"  or  "extended" 
integral  equations. 


H.  (r) 

inc 


j  K(r')  x  V'v|/(  r  ,r  '  )dS  ' 
S 


(r  inside  S) 


E.  (  r )  = 
inc 


1 

Tore 


k2  Ky  -  (V'  K)V> 


dS', 


(la) 


(lb) 


where  k(=to/c=2nA)  is  the  free-space  propagation  constant,  £Qis  the  permittivity  of  free 

space,  and  \\i  (r,r')  is  the  free-space  Green's  function  exp  (ik |  r-r' |  )/4rc|  r-r'  ] . 

Let  the  observation  point  r  in  Eqs.  (la)  and  (lb)  approach  the  surface  S  of  the 
conductor  from  inside  S,  and  convert  the  surface  integrations  in  Eq.  (1)  to  circular 
principal-value  integrations  using  the  following  formula  (derived  by  a  straightforward 

g 

integration  near  the  singularity  of  \ \i  ): 


V>  dS'  =  |  V>  dS'  - 
S 


A 

n 


S(r  >S  ) 


(2) 


where  o  denotes  the  principal- value  surface  integration  evaluated  by  excluding  the 


S 

singular  point,  r'=  r.  of  the  integrand  by  a  limiting  circular  "principal  area"  centered 

on  r,  and  n  is  the  outward  unit  normal  from  the  surface  S  at  r.  Eqs.  (1)  then 
yield  the  augmented  magnetic  arm  electric-field  surface  integral  equations  for  the 
exterior  scattering  problem: 
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(3a) 


-  Hinc(r)  =  -i-  n  x  K  +  j  K  x  V>  dS' 

S 


( r  on  S ) 


E  (r)  = 
me 


(V 


s 

2  1(0  £ 


K)n 


I0)£ 


k2  K y  -  (V'  K)V>jdS'. 


(3b) 


Taking  ft  cross  these  overdetermined  Eqs.  (3)  reduces  them  to  the  original  even 
determined  MFIE  and  EFIE, 


ft  x  H. 


inc 


y  K  -  ft  X  j  K  x  V>  dS' 
S 


(r  on  S) 


ft  x  E. 


inc 


A 

n 

1C0£ 


k2  Kv|/ 


(V'K)V> 


dS' 


(4a) 


(4b) 


As  mentioned  in  the  Introduction,  the  original  integral  Eqs.  (4)  are  plagued  by 

spurious  solutions  for  K  at  the  resonant  frequencies  of  the  cavity  formed  by  the  surface 
S.  Although  either  of  the  augmented  Eqs.  (3)  eliminate  the  spurious  solutions  for  most 
shapes,  they  must  both  be  used,  in  general,  to  eliminate  all  the  spurious  resonances 
when  the  surface  S  is  a  body  of  revolution.  The  combined-field  integral  equation9  10  11 
eliminates  the  spurious  resonances  by  adding  Eq.  (4a)  to  -a^ft  crossed  into  Eq.  (4b),  and 
uniqueness  of  solution  of  the  combined-source  equation  “  1  "  14  follows  from  its  operator 
being  the  adjoint  of  the  combined-field  operator"'  14  (The  real  constant  is  often 
chosen  equal  to  the  free-space  wave  admittance.) 

To  derive  the  dual-surface  integral  equations,  return  to  the  extended  integral 

equations  (1).  The  current  K(r)  in  Eq.  (la)  or  (lb)  is  uniquely  determined  at  every 
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frequency  if  Eq.  (la)  or  (lb)  is  satisfied  for  all  r  inside  S16.  Conceivably,  one  could 

determine  the  current  K  by  solving  numerically  the  vastly  overdetermined  set  of  extended 

integral  equations  that  results  from  Eq.  (la)  or  (lb)  applied  to  points  r  separated  by  a 
small  fraction  of  a  wavelength  throughout  the  volume  enclosed  bv  S.  Or  one  could 
supplement  the  surface  integral  equations  (4)  with  the  extended  integral  equations  (1) 

applied  at  selected  points  r  within  S.l7,is'w  The  former  approach  introduces  a 
prohibitive  number  of  equations  for  multi-wavelength  bodies,  and  in  applying  the  latter 
approach  one  has  no  convenient,  reliable  criterion  for  selecting  the  number  and  position 
of  the  interior  points  at  which  the  extended  integral  equations  (la)  and  (Ib)must  be 

satisfied  to  assure  Eqs.  (4a)  and  (4b)  produce  the  correct  unique  current  K  at  all 
frequencies.  (The  modified  Green’s  function  method"0'*1'""  for  eliminating  the  spurious 
solutions  from  the  original  surface  integral  equations  (4)  suffers  from  a  similar 
uncertainty  in  choosing  the  proper  number  and  origin  of  eigenfunctions  in  the 
representation  of  the  modified  Green  s  function."  ) 


|h  Waterman,  P  C.  (1965)  Matrix  formulation  of  electromagnetic  scattering.  Proc.  IEEE , 
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If,  however,  the  extended  integral  equations  (la)  and  (lb)  are  incorporated  at  points 

r  on  a  surface  Sg  parallel  to,  and  a  small  distance  8  >  0  inside  the  current  surface  S 
of  the  perfect  conductor  (see  Figure  1),  by  adding  an  cross  Eqs.  (la)  and  (lb)  at  these 
points  to  the  original  MFIE  Eq.  (4a)  and  EFIE  Eq.  (4b),  respectively,  one  obtains  the 
"dual-surface"  magnetic  and  electric-field  integral  equations: 

(5a) 

n  x  HQ(r)  =  y-  K(r)  -  ft  x  o  K(r')  x  V'\j/o(r,r')dS' 


n  x  E  (r)  =  A  x  f  Ik2  Ky  -  (V'  K)V>  IdS', 
ov  itoe  To  s  To 

O  o 


where  E  ,  H  and  ur  are  defined  as 
o  o  To 


E  (r)  s  E  (r)  +  a  E.  (r  -  8n) 
ov  tncv  mcv 


H  (r)  =  H.  (r)  +  a  H-  (r  -  5n) 
o  mcv  tncv 


(r  on  S) 


V|fo(r,r')  =  V|/(r,r')  +  a  vj/(r  -  8ft,  r'). 
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Figure  1.  Geometry  of  a  Perfect  C  iuctor  with  Current 
Surface  S  and  Parallel  surface  Sg. 


These  dual-surface  magnetic  and  electric-field  integral  equations,  (5a)  and  (5b),  although 
identical  in  form  and  comparable  in  complexity  to  the  original  MF1E  (4a)  and  EFIE  (4b), 

provide  a  unique  solution  for  K  at  all  real  frequencies  as  long  as  the  constant  a  is 
imaginary  and  the  positive  real  constant  5  is  less  than  about  A/2.  (In  the  numerical 
solutions  of  Eq.  (5a)  described  in  Section  4  below  we  choose  a  equal  to  i  and  5  equal  to 
the  smaller  of  about  A/4  or  1/4  the  breadth  of  the  scatterer  along  the  normal  at  the 

point  r.  Using  a’s  of  ±.5i,  ±i,  and  ±1.5i.  and  varying  5  from  A/8  to  3A/8  did  not 
significantly  change  the  computed  solution,  although  the  number  of  iterations  required 
by  the  conjugate  gradient  method  to  attain  the  same  value  of  normalized  residual  error 
varied  somewhat  with  a  and  8.) 
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3.  UNIQUENESS  OF  SOLUTION  OF  THE  DUAL-SURFACE  INTEGRAL 
EQUATIONS 

Uniqueness  of  solution  for  the  dual-surface  magnetic-  and  electric-field  integral 
equations  can  be  proven  by  considering  the  fields  radiated  by  the  solution  currents. 

Concentrating  on  the  dual-surface  MFIE  first,  let  Hs(r)  be  the  magnetic  field  radiated 

by  the  solution  K  to  Eq.  (5a);  specifically 


Hs(r)  =  J  K(r)  x  V'\y(r,r')dS'  (r  not  in  S). 
S 


(7) 


If  K  were  the  correct  unique  current  for  this  scattering  problem,  Hs(r)  in  Eq.  (7)  would 

be  the  correct  scattered  field  for  all  r  not  in  the  surface  current.  However,  since  we 

do  not  know  at  this  point  that  the  solution  K  to  Eq.  (5a)  is  the  correct  unique 

solution,  Eq.  (7)  simply  defines  an  unknown  magnetic  field  Hs(r). 

Taking  the  curl  of  Eq.  (7)  twice  reveals  that  this  unknown  magnetic  field  satisfies 

the  homogeneous  vector  wave  equation  for  all  r  not  in  the  surface  current  K,  that  is, 

V  x  V  x  H  -  k2Hs  =  0  (r  not  in  S).  (8) 

Letting  r  approach  S  in  Eq.  (7)  from  the  inside  of  S,  and  using  the  principal -value 

formula  [Eq.  (2)J,  we  obtain 


Hs(r-)  =  x  K(r)  +  <j>  K(r')xV'tj/(r,r')  dS'  (r  on  S),  (9) 

S 

where  r-  in  Hs(r-)  indicates  the  field  evaluated  just  inside  the  surface  current.  Since 
Eq.  (7)  holds  for  all  r  inside  S,  we  can  express  Hs  on  the  parallel  surface  Sg  as 
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Add  Eq.  (9)  to  Eq.  (It))  multiplied  by  a  and  take  n  cross  the  result,  to  get 


ft  x 


(H  (r-)  +  uH  (r  -  5n) 

s  s 


K(r)  +  ft  x  j,  K(r')  x  V>q(  r,r')dS'.  (II) 
^  (r  onS) 


Comparing  Eq.  (11)  with  Eq.  (5a),  which  K  must  aJso  satisfy,  reveals 


ft  x 


(Hs(r-)  +  Hjnr( r ))  +  a  (Hs(r-5n)  +  Hjnc(r-6n)) 


=  0  (r  on  S). 


(12) 


The  incident  magnetic  field  also  satisfies  the  vector  wave  equation 

V  x  V  x  H.  -  k:  H.  =  0  .  (13) 

me  me 


Thus,  we  can  add  Eq.  (8/  to  Eq.  (13).  and  rewrite  Eq.  (12)  to  arrive  at  the  interior 
boundary  value  problem, 

V  x  V  x  H  -  k~  H  =  0  (r  inside  S) 


A 

n 


H(r)  +  a  H(r 


5ft) 


=  0 


(r  -*S  from  inside), 


(14a) 

(14b) 


where  the  total  magnetic  field  H(r)  is  given  by  the  sum  of  H^,c(r)  and  Hs(r). 

The  final  steps  of  the  uniqueness  proof  consist  in  showing  that  the  boundary  value 

problem  defined  by  Eq.  (14)  has  only  the  trivial  solution,  H(r)  =  0,  for  the  total  field 
throughout  the  volume  enclosed  by  the  surface  S,  provided  the  constant  cl  is  imaginary 


and  the  positive  real  constant  6  is  smaller  than  about  X/2.  To  show  this,  rewrite  the 
boundary  condition  of  Eq.  (14b)  explicitly  for  the  magnetic  field  tangent  to  the 
surface  S 


H{(r)  +  a  H{(r  -  8n)  =  0  (r-^S  from  inside).  (15) 

The  tangential  magnetic  fields,  H((r)  ami  H{(r  -  8n),  are  complex  numbers  that  can  be 
expressed  in  the  form  of  a  magnitude  and  phase 

(16a) 

Ht(7)  =  |Ht(7)|ei<Kr) 

(r-^S  from  inside) 

H{(7  -  8n)  =  f  |  Ht(7) }  +  AHt(7,5)Je1t<t>(r)  +  A<t)(r*5)]  (16b) 


where  AHf  and  A<)>  are  the  differences  between  the  magnitudes  and  phases  of  Ht(r)  and 
H((r  -  6n).  Insert  Ht(r)  and  Ht(r  -  5n)  from  Eq.  (16)  into  Eq.  (15)  to  get 


+  AH 


cos  A<}>  +  i  sin  A<|> 


=  0. 


(17) 


Because  |  H{  | ,  AHt,  and  A<J)  are  real  numbers,  if  we  let  the  constant  a  be  an  imaginary 
number(iaj),  the  real  and  imaginary  parts  of  Eq.  (17)  equate  separately  to  give 


a. 

i 


|  +  AH(  sin  A<j>  =  0 


(18a) 


+  AH(  cos  A<j>  =  0 


(18b) 
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o  + 

For  small  0,  A<(>  will  be  small  -  certainly  not  ±  90  -  and  thus  Eqs.  (18)  imply 

Hf  =  0  and  AH(  =  0,  or 

H  (r)  =  0  (19a) 

(r  *S  from  inside) 

Ht(r  -  5n)  =  0.  (19b) 

In  other  words,  when  the  constant  a  is  chosen  imaginary  and  5  is  not  large,  the  two 
separated  tangential  fields  in  the  boundary  condition  of  Eq.  (14b)  are  each  zero.  That 
is,  the  tangential  magnetic  field  on  both  the  surface  S  and  Sg  are  zero. 

The  boundary  condition  of  Eq.  (19a)  restricts  the  nonzero  solutions  of  Eq.  (14a)  to 
the  resonant  modes  of  the  cavity  formed  by  a  perfectly  magnetically  conducting  surface 
S.  These  modes,  which  exist  for  a  given  cavity  only  at  discrete  frequencies,  fonn 
standing  waves  within  the  cavity  with  magnetic  and  electric  fields  that  can  be  chosen 
retil  and  imaginary,  respectively.'4  In  particular,  the  tangential  magnetic  field  near 
the  surface  S  can  be  expressed  approximately  as 

H  ( r  .  r  )  ~  A(r  ,  r  )  sin  yr  (20a) 

t  s  n  s  n  *  n 

where  (r  .rn)  are  the  coordinates  tangent  and  normal  to  the  surface  S,  y  is  a  positive 
real  propagation  constant  with  a  value  equal  to  or  less  than  the  propagation  constant  k 

of  free  space,  and  the  amplitude  A(rs>  r  )  varies  with  r^slowly  compared  to  the 
variation  of  sin  yr  (With  respect  to  the  r^  direction  the  cavity  can  be  considered  a 
shorted  waveguide  with  varying  cross  section.)  If  we  let  r  =  0  on  the  surface  S,  the 


*The  one  exception  would  be  if  there  were  a  zero  of  H(  near  the  surface  S  or  Sg.  In 

that  case  we  can  expand  the  boundary  condition  of  Eq.  15  along  the  normal  direction  r^ 

in  a  Taylor  series  about  the  zero  to  show  that  flH^dr  must  also  vanish  at  the  zero  of  H{ 

for  imaginary  a  and  small  5.  Since  Eq.  (20a)  shows  that  no  cavity  can  support  modes 
with  both  die  tangential  field  and  its  normal  derivative  zero  on  its  surface,  the 
solution  to  Eq  (14)  is  unique  in  this  exceptional  case  as  well. 

24  Borgnis,  F.E  and  Papas,  C.H.  (1958)  Electromagnetic  waveguides  and  resonators,  E/tcvclo- 
pi’dia  of  Physics ,  16  Ed.  S.  Flugge,  Springer-Verlag,  Berlin. 
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boundary  condition  of  Eq.  (19a)  is  satisfied  by  (20a).  The  boundary  condition  of 
Eq.  (19b)  applied  to  (20a)  requires  that 


y5  =  m  7i  (20b) 

for  m  equal  to  a  positive  integer.  (We  assume  that  there  will  be  some  portion  of  the 
surface  S  where  A  will  not  be  zero.  For  if  the  tangential  magnetic  field  were  zero 

throughout  the  volume  between  S  and  Sg,  the  fields  would  be  zero  throughout  the  cavity. 

Also,  if  degenerate  modes  exist,  we  assume  their  H(  fields  will  be  linearly  independent 
over  the  surface  Sg  and  thus  Eq.  (20b)  will  still  hold.  Because  the  maximum  value  of  y 

is  k  =  2nA,  the  condition  of  Eq.  (20b)  cannot  be  satisfied  for 

0  <  6  <  A/2  .  (21) 

The  approximate  sign  is  included  in  the  right  side  of  the  inequality  in  (21)  because 
(20a)  is  an  approximate  expression  for  the  standing  wave  field  near  the  surface.  If  we 
look  specifically  at  the  resonant  cavity  formed  by  shorting  the  ends  of  a  waveguide  of 
arbitrary  uniform  cross  section,  we  find  (20a)  applies  exactly  with  A  equal  to  a 
constant.  Thus,  the  inequality  of  (21)  holds  exactly  for  a  shorted  waveguide  cavity. 

For  a  spherical  cavity  the  fields  vary  radially  as  spherical  Bessel  functions  of  the 
first  kind.  For  asymptotically  large  spheres,  (20a)  again  holds  exactly,  and  for  all 
spheres  large  enough  to  sustain  resonant  modes,  (20a)  holds  to  a  good  approximation  near 
the  surface  -  thereby  confirming  the  approximate  inequality  [Eq.  (21)]  for  spherical 
cavities. 

In  summary,  the  only  solution  to  Eq.  (14)  for  a  imaginary  and  0  <  5  <  A/2  is  the 

trivial  solution,  H(r)  =  H^c(r)  +  H(r)  =  0  throughout  the  volume  enclosed  by  S.  Since 

E  =-V  x  H/itoe  ,  the  electric  field  E(r)  =  E.  (r)  +  E(r)  within  this  volume  is  also 
o  me  « 

identically  zero.  And,  as  mentioned  in  Section  2,  it  is  a  simple  matter  to  prove  that 
the  current  that  produces  the  negative  of  the  incident  electromagnetic  fields  throughout 
the  volume  enclosed  by  S  is  the  correct  unique  current  for  the  exterior  scattering 

problem.  (Namely,  E  and  H  equaling  zero  inside  S  implies  n  x  E  =  0  and  n  x  H  =  K  for 
the  fields  just  outside  S  —  the  conditions  required  for  uniqueness  of  solution  of  the 
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exterior  problem'5. )  Since  this  unique  solution  has  been  derived  from  the  solution 
current  of  Fq  (5ai.  the  dual-surface  magnetic-field  integral  equation  (5a)  has  a  unique 
solution. 

Beginning  with  the  solution  current  of  the  dual-surface  electric-field  integral 
equation  (5b).  and  defining  the  electric  field 


Em 


s 


I 

i  roc 

o 


j  k:  Ky 

S 


<V'- 


K)Vy 


dS' 


( r  not  in  S ) 


(22) 


initially,  instead  of  the  magnetic  field  Eq.  (7).  we  obtained  the  same  inequality  (21)  as 
the  sufficient  condition  for  the  uniqueness  of  solution  of  the  dual-surface 
electric-field  integral  equation  (5b). 

In  numerical  practice,  we  suggest  choosing  «  equal  to  i  to  weight  the  fields  on  S  and 
S -  equally  in  the  boundary  condition  |Eq.  ( 1 4 b ) J  by  an  imaginary  constant.  Likewise,  we 
suggest  choosing  5  equal  to  about  X/4,  to  keep  the  surface  Sg  about  tut  equal  distance 
between  the  two  critical  values.  5  =  0  and  X/2  The  dual-surface  integral  equations 
allow  spurious  solutions  at  5  =  0  where  they  reduce  to  the  original  integral  equations 
(4),  and  at  5  equal  to  or  greater  than  about  X/2  where  the  dual -surface  boundary 
condition  | Hq.  (14b)]  no  longer  insures  uniqueness  of  solution.  When  the  breadth  of  the 
scatterer  along  the  normal  is  less  than  X.  one  can  choose  5  equal  to  1/4  the  breadth 
instead  of  X/4. 

A  numerical  demonstration  of  the  elimination  of  the  spurious  resonances  by  the 
dual  surface  magnetic- field  integral  equation  is  given  in  Figures  2  and  3.  Figure  2 
plots  the  total  (integrated)  radar  cross  section  versus  the  perimeter  of  a  perfectly 
conducting  cube  of  side  length  s  as  computed  using  the  conventional  magnetic  field 
integral  equation  (4a).  The  spurious  resonances  begin  to  contaminate  the  MFIE  solution 
m  Figure  2  near  the  first  resonance  of  the  cube  at  4s/X  =  2.8  and  continue  to  distori 
the  solution  at  an  increasing  rate  commensurate  with  the  increasing  density  of  resonant 
frequencies.  Figure  3  shows  clearly  that  the  dual-surface  magnetic-field  integral 
equation  (5a)  eliminates  the  spurious  resonances  from  the  MFIE  solution  in  Figure  2. 


Muller,  C  (1969)  Foundations  of  the  Mathematical  Theory  of  Electromagnetic  Waves. 
Springer-Verlag.  New  Yo:!;,  Theorem  71 
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Figure  2.  Total  Radar  Cross  Section  Versus  Perimeter  of  a 
Perfectly  Conducting  Cube  as  Computed  with  the  Conventional 
Magnetic-Field  Integral  Equation  (4a). 


16 


M 


Figure  3.  Total  Radar  Cross  Section  Versus  Perimeter  of  a 
Perfectly  Conducting  Cube  as  Computed  with  the  Dual-Surface 
Magnetic-Field  Integral  Equation  (5a). 
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4.  NUMERICAL  SOLUTION  TO  THE  DUAL-SURFACE  MAGNETIC- 
FIELD  INTEGRAL  EQUATION  BY  THE  CONJUGATE 
GRADIENT  METHOD 

The  similarity  of  the  dual-surface  integral  equations  (5)  to  the  original  integral 
equations  (4)  allows  them  to  be  solved  numerically  by  a  minor  modification  to  existing 
MFIE  and  EFIE  computer  programs.  One  merely  adds  the  values  of  the  incident  field  and 

free-space  Green’s  function  (each  multiplied  by  a)  at  the  points  r  -  5n  to  their 

respective  values  at  r  used  in  the  computer  programs  of  the  original  integral  equations. 

In  particular,  we  consider  ?  straightforward  numerical  solution  to  the  dual-surface 
magnetic-field  integral  equation  (5a)  for  scattering  from  the  perfect  conductor  S. 

Divide  the  surface  S  of  the  scatterer  into  M  patches,  assume  the  current  is  a 
constant  vector  over  each  patch,  approximate  the  value  of  the  Green’s  function  V'y  over 
each  patch  by  a  constant  vector  equal  to  the  value  of  V'\j /  at  the  center  of  the  patch, 
and  apply  the  integral  equation  (5a)  at  the  center  of  each  patch.  In  short,  approximate 
the  integral  in  Eq.  (5a)  by  the  summation 


n-  x  H  (r) 

1  o  i 


where  AS.  is  the  area  of  each  patch,  and  the  self-patch  (i  =  j)  in  the  summation  is 

J  t  g 

taken  as  the  "principal  area”  excluded  by  the  principal-value  integral  in  Eq.  (5a).  (In 
the  language  of  the  method  of  moments,  we  have  used  puls^  basis  functions  and  delta 
testing  functions.) 

For  each  patch  there  are  two  complex  unknown  components  of  surface  current  K  and 
two  complex  scalar  equations.  Thus,  Eq.  (23)  represents  a  simultaneous  set  of  2M 
linear  complex  equations  for  2M  complex  unknowns,  and  can  be  written  in  tensor  notation 
as 


a.,  x. 


«J  J 


i  =  1,2 . N  =  2M, 


(24) 
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where  the  Xj  are  the  complex  unknown  components  of  surface  current,  the  a^.  are  the 
elements  of  the  given  coefficient  matrix,  and  the  b.  are  the  given  incident-field 
values.  Summation  from  1  to  N  over  the  repeated  index  j  in  Eq.  (24)  is,  of  course, 
implied. 

In  solving  Eq.  (24)  for  three-dimensional  bodies  on  readily  available  computers,  one 
quickly  encounters  the  problem  of  limited  central  or  virtual  memory  and  excessive 
computer  time  as  the  size  of  the  body  is  increased  beyond  a  wavelength.  For  example, 
we  have  found  that  the  dual-surface  MFIE  requires  a  minimum  of  about  25  patches  per 
square  wavelength*  to  achieve  reasonable  accuracy  in  the  computed  currents  and  far 
fields.  A  cube  of  side  length  s  thus  requires  about  M  =  150  (s fk)~  patches  or 


N  =  2M  =  300(sA): 


(25) 


complex  unknowns,  and  a  memory  of 


W  =  2N2 


(26) 


real  words  W  to  store  the  complex  N  x  N  matrix  with  elements  a^.  On  our  VAX  8650 
with  an  allotted  virtual  memory  of  one  million  words,  the  solution  of  Eq.  (24)  is 
limited  to  cubes  less  than  about  1.5  wavelengths  on  a  side.  (Paging  time  became 
excessive  as  the  amount  of  virtual  memory  used  approached  a  million  words;  and  allotting 
more  virtual  memory  to  the  computer  was  not  as  efficient  an  alternative  as  using  direct 
access  files  for  increasing  computer  storage  capacity.'6)  Moreover,  we  found  that 


*  This  requirement  of  about  5  linear  divisions  per  wavelength  to  get  reasonable 
accuracy  is  not  surprising  if  one  considers  that  the  current  and  the  Green’s  function  in 
the  surface  integral  equation  varies  along  the  surface  of  the  scatterer  with  a  maximum 
spatial  frequency  equal  to  about  one  cycle  per  wavelength.  This  means  that  the  product 
of  the  current  and  Green’s  function  in  the  integral  of  the  integral  equation  has  a 
maximum  spatial  frequency  of  about  2  cycles  per  wavelength.  The  sampling  theorem 
would  then  require  about  4  samples  per  linear  wavelength  to  accurately  approximate  the 
integral  of  the  current  times  the  Green's  function  by  a  summation. 


Perry,  T.S.  and  Zorpette,  G.  (1989)  Supercomputer  experts  predict  expansive  growth, 
IF.FJi  Spectrum ,  26:  26-33. 
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solving  Eq.  (24)  using  Gaussian  elimination  on  this  computer  with  a  32-bit  Linpack 
performance  rating5  of  1.3  megaflops  took  a  CPU  time  Y  given  approximately  (for  large 
N)  by 

T  =  2.2  x  10 1  N3  (27) 

GE 

minutes,  which  becomes  prohibitive  for  cubes  greater  than  about  2  wavelengths  on  a  side. 
Incorporating  two-fold  symmetry  of  the  cube  to  reduce  N  by  a  factor  of  4  still  did  not 
allow  us,  in  a  reasonable  computer  time,  to  solve  for  scattering  from  cubes  larger  than 
about  4  wavelengths  on  a  side  using  Gaussian  elimination  on  the  available  computer.  The 
fastest  available  computers  (those  with  32-bit  Linpack  performance  ratings  of  about  100 
megaflops)  would  take  many  hours  of  computer  time  to  solve  Eq.  (24)  using  Gaussian 
elimination  for  general  3-D  scatterers  larger  than  a  few  wavelengths  across.  Even  with 
massive  vector  parallelism  it  is  difficult  to  conceive  of  digital  computers  extending 
appreciably  the  formidable  restriction  on  sA  presented  by  the  (sA)6  dependence  of  the 
computer  time  in  Eq.  (27)  for  solving  Eq.  (24)  using  Gaussian  elimination. 

To  extend  the  limits  of  computer  storage  and  processing  time  on  the  available 
mainframe,  we  made  use  of  direct  access  memory  files  on  disk  and  solved  Eq.  (24) 
iteratively  using  the  conjugate  gradient  method  rather  than  Gaussian  elimination.1  A 
direct  access  file  was  used  to  store  the  rows  of  the  N  x  N  complex  matrix.  The  file 
could  be  either  opened,  written  to,  or  read  from  by  a  single  Fortran  command,  and 
increased  our  available  memory  from  1  million  to  30  million  words.  (Of  course,  with 
iterative  solvers  one  can  greatly  reduce  computer  storage  requirements  by  generating  the 
coefficient  matrix  during  each  iteration.  However,  this  greatly  increases  the  required 
CPU  time.) 

The  drawbacks  of  using  direct  access  files  are  the  necessary  additional  computer 
programming,  the  somewhat  greater  CPU  time,  and  possibly  a  large  increase  in 
input/output  time.  Use  of  direct  access  files  roughly  doubles  the  CPU  time  required  to 
solve  Eq.  (24)  using  the  conjugate  gradient  method.  Gaussian  elimination  CPU  tunes  are 
either  roughly  doubled  or  multiplied  by  a  factor  of  about  10,  when  using  direct  access 
files,  depending  on  whether  or  not  a  round-off  error  check  is  included  in  the  Gaussian 
elimination  algorithm.  (We  shall  discuss  this  later  in  conjunction  with  Tables  1 
and  2.)  The  extra  input/output  time  associated  with  the  direct  access  files  may 
dominate  computer  turn-around  time  on  our  computer  system  when  the  matrix  is  solved 
using  the  conjugate  gradient  method.  Meaningful  input/output  times  are  elusive, 
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however,  on  central  computers  since  they  depend  so  strongly  on  the  particular  direct 
access  system,  the  way  the  system  is  tuned,  turd  the  number  of  users  sharing  the  machine 

We  found  the  conjugate  gradient  iterative  method  tut  efficient  alternative  to  Gaussitui 
elimination  for  solving  the  large  system  of  equations  (24)  generated  by  the  dual-surface 
magnetic-field  integral  equation  (5a)  applied  to  3-D  multi-wavelength  bodies.  It 
converges  in  a  finite  number  of  steps  for  tiny  initial  guess,  as  long  as  the  matrix  is 
not  singular  and  round-off  errors  are  kept  negligible.  The  flexibility  of  accepting  any 
initial  guess  permits  the  user  to  stop  the  iteration  and  start  using  the  current 
estimate  for  the  solution  as  the  new  initial  guess.  This  restarting  eliminates  the 
accumulated  round-off  errors,  but  requires  ;m  extra  matrix  multiplication,  slows 
convergence,  and  increases  the  run  time  Thus,  resuming  should  be  done  only  when 
necessary'  to  reduce  the  round-off  errors.  For  solving  the  dual-surface  equations 
(23-24)  on  our  32-bit  computer,  it  was  sufficient  to  restart  only  once  before  the  final 
step,  even  when  N  was  as  large  as  3468.  (Wang  and  Dubberley"  discuss  this  restart 
problem. ) 

We  used  a  version  of  the  conjugate  gradient  method  referred  to  as  "Case  A"  in 
Reference  2.  We  applied  the  conjugate  gradient  algorithm  to  the  matrix  with  elements 
£i- 1  in  Eiq.  (24)  rather  than  directly  to  the  operator  of  ,.e  uimI-s. inace  magnetic-field 
integral  equation  (5a).  Some  justification  this  may  he  found  in  the  recent 
papers'*  where  it  is  concluded  that  iterative  techniques  applied  directly  to  the 
operator  and  implemented  numerical. y  con...!,  c  imp'icb  discretization  and.  in  many 
cases,  a  corresponding  moment  method  interpretation.  It  is  interesting  that,  as  early 
as  1443.  Hotelling  "  in  his  review  of  some  new  methods  in  matrix  calculation  commented 
that,  "The  combination  of  this  device  |Mallock  electrical  calculating  machine]  with  the 
iterative  method... 


W.mg.  J  J  11.  and  Dubberley,  J.R  (14,34)  Computation  of  electromagnetic  fields  in  large 
biological  bodies  by  an  iterative  moment  method  with  a  restart  technique,  IEEE  Trans. 
Microwave  Theory  and  Techniques.  MTT-37:  1418-1423 

Ray,  SI-  and  Peterson.  A  F.  (1488)  Error  and  convergence  in  numerical 
implementations  of  the  conjugate  gradient  method.  IEEE  Trans  Antennas  Proposal.  36: 
1824-1827. 

Sarkar,  I  K  .  Yang.  X  .  and  Areas.  F.,  (1488)  A  limited  survey  of  various  conjugate 
gradient  methods  for  solving  complex  matrix  equations  arising  in  electromagnetic  wave 
interactions.  Wave  Motion  10(6):  527-546 

Hotelling.  H  (1443)  Some  new  methods  in  matrix  calculation.  Annals  of  Math.  Stat  . 
I4(h:  1-34. 
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offers  what  seems  at  present  [1943]  the  best  hope  for  the  systematic  inversion 
[solution]  of  large  matrices.'1 

Table  1  compares  the  number  of  major  complex  operations  required  to  solve  large 
matrices  by  means  of  Gaussian  elimination  and  the  conjugate  gradient  method,  when  using 
direct  access  memory  files.  Table  2  shows  the  associated  CPU  times  required  by  our 
32-bit,  1.3  megaflop  Linpack-performance-rated  VAX  8650  computer. 


iable  1.  Number  of  Complex  Operations  Required  for  N  x  N 
Matrix  Solution  Using  Gaussian  Elimination  and  the  Conjugate 
Gradient  Method. 


METHOD 

NUMBER  OF  OPERATIONS 

| 

Add 

Subtract 

Do  Loops 

Elements 

Written 

Elements 

Read 

If 

2N3 

0 

In3 

In3 

In3 

3-N3 

In3 

3 

3 

3 

2 

2 

3 

Conjugate  Grad, 
(per  iteration) 

2N2 

2N2 

0 

2N2 

0 

2N2 

0 

Table  2.  CPU  Time  Required  for  Complex  Operations  on  Our  32-Bit, 
1.3  Megaflop  Linpack-Performance-Rated  Computer. 


OPERATION 

CPU  Time  (10'®  sec) 

Complex  Add 

0.84 

Complex  Subtract 

0.96 

Complex  Multiply 

2.03 

Complex  Divide 

13.09 

if 

15.09 

Do  Loop 

0.67 

Read  per  complex  element  (for 
large  N) 

2.30 

Write  per  complex  element  (for 
large  N) 

2.41 
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The  CPU  times  for  the  complex  operations  in  Table  2  can  be  inserted  into  Table  1  to 
estimate  the  total  CPU  times,  t  amt  t  ,  for  this  computer  to  solve  Eq.  24  by  Gaussian 
elimination  and  the  conjugate  gradient  method.  Specifically, 


t  =  1.93  x  10 7  N1  (28) 

( if: 

t  =  1.95  x  IQ'7  N:I  (29) 


minutes,  where  N  is,  as  usual,  the  dimension  of  the  complex  matrix,  and  I  is  the  number 
of  iterations  needed  for  convergence  using  the  conjugate  gradient  method.  Comparing  the 
estimated  CPU  time  [Eq.  (28) j,  for  the  solution  to  scattering  from  the  cube  by  Gaussian 
elimination,  with  the  actual  CPU  run  times  (for  the  whole  program)  given  approximately 
by  the  formula  (27),  one  finds  that  the  estimated  time  [Eq.  (28)]  is  about  90  percent  of 
the  actual  total  CPU  run  times.  Likewise,  Eq.  (29)  gives  an  estimated  CPU  time  for  the 
conjugate  gradient  method  that  is  about  75  percent  of  the  actual  total  CPU  times  for 
scattering  from  large  cubes  (see  Table  3).  The  additional  10  percent  and  25  percent  CPU 
tunes  are  taken  mainly  by  matrix-fill,  complex  conjugate,  and  miscellaneous  overhead 
operations. 


Table  3.  Number  of  Iterations  (I)  and  Actual  Total  CPU  Time  Using  the 
Conjugate  Gradient  Method  on  our  32-Bit,  1.3  Megaflop  Linpack- 
Performance- Rated  Computer. 


s  A 

N 

PatchesA2 

I 

1/N 

CPU  Time 
(h:m:s) 

0.75 

■ 

28 

35 

.73 

0:00:04 

0.75 

63 

39 

.36 

0:00:12 

0.75 

1 12 

42 

.22 

0:00:32 

■ 

28 

61 

.32 

0:00:45 

B 

1 

63 

62 

.14 

0:03:08 

1 13 

61 

.08 

0:10:47 

WMSMm 

■ 

25 

83 

.19 

0:04:29 

44 

82 

.11 

0:13:53 

1 

69 

88 

.07 

0:35:15 

768 

28 

90 

.12 

0:13:30 

1200 

44 

92 

.08 

0:33:21 

1728 

63 

93 

.05 

1 :09:44 

5.0 

1728 

23 

118 

.07 

1:25:47 

5.0 

3468 

46 

119 

.03 

6:23:10 

6.75 

3468 

25 

141 

.04 

7:28:34 
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The  time  estimates  Eqs.  (28)  and  (29)  reveal  that,  for  our  typical  computer,  the 
solution  to  Eq.  (24)  by  the  conjugate  gradient  method  will  take  less  computer  time  than 
Gaussian  elimination  if  the  number  of  required  iterations  I  is  less  than  N.  This 
conclusion  holds  whether  or  not  direct  access  files  are  used,  because  including  the 
write  and  read  statements  from  Table  1  due  to  the  use  of  direct  access  files  roughly 
doubles,  as  mentioned  above,  the  computer  times  for  both  Gaussian  elimination  (with 
round-off  error  check)  and  the  conjugate  gradient  method. 

It  is  important  to  note,  however,  that  the  logical  IF  operations  and  one  half  the 
multiplication  operations  listed  in  Table  1  for  Gaussian  elimination  are  produced  by  the 
round-off  error  check  in  our  Gaussian  elimination  algorithm.  If  this  round-off  error 
check  is  omitted,  the  revised  Gaussian  elimination  CPU  time  t'  estimated  from  Tables  1 

GE 

and  2  is  given  by 


t;E  =  10-7  N3  (30) 

minutes,  about  one  half  the  CPU  times  given  by  Eq.  (27)  or  (28).  Comparing  Eq.  (30) 
with  Eq.  (29)  shows  that  the  conjugate  gradient  method  becomes  faster  than  Gaussian 
elimination  without  the  round-off  error  check  when  the  number  of  iterations  I  is  less 
than  about  N/2. 

If  in  addition  to  omitting  the  round-off  error  check  from  the  Gaussian  elimination 
algorithm,  all  computations  could  be  done  in  central  memory  without  using  direct  access 
files,  the  write  and  read  operations  would  be  eliminated  from  Table  l,  and  CPU  times 
[Eqs.  (28)  and  (29)1,  would  be  replaced  by 


t°  =  0.2  x  10  7  N3 

GE 


(31) 


t°  =  1.2  x  10 7  N2! 


(32) 


minutes.  Comparison  of  Eqs.  (31)  and  (32)  reveals  that,  if  all  computations  can  be 
handled  in  central  memory,  the  conjugate  gradient  method  is  faster  than  Gaussian 
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elimination  (without  a  round-off  error  check)  when  the  number  of  required  iterations  is 
less  than  about  N/6.'M 

Table  3  lists  the  number  of  iterations  and  actual  CPU  tunes  using  the  conjugate 
gradient  method  for  the  dual-surface  magnetic-field  equations,  (23-24),  applied  to 
plane-wave  scattering  from  a  perfectly  conducting  cube.4  The  plane  wave  was  incident 
broadside  upon  the  cube  of  side  length  s,  and  the  parameters  a  and  5  in  the  dual-surface 
magnetic-field  integral  equation  (5a)  were  set  equal  to  i  and  '/if>  A.,  respectively.  The 
initial  value  taken  for  the  solution  vector  in  the  conjugate  gradient  algorithm  was 
zero,  and  the  iterations  were  terminated  when  the  ratio  of  the  magnitude  of  the  residual 
vector  to  the  magnitude  of  the  source  vector  became  less  than  10  *.  Two-fold  (xy) 
symmetry  of  the  cube  was  used  to  reduce  the  number  of  unknowns  N  in  Table  3  by  a  factor 
of  4,  so  that  N  equals  75(s AT  rather  than  300(sA)*  [see  Eq.  (25)]  for  N  in  terms  of 
the  side-to-wavelength  ratio  (s A)  of  a  cube  with  25  patches  per  square  wavelength,  the 
minimum  number  needed  to  achieve  reasonable  accuracy  in  the  computed  surface  currents 
and  far  fields. 

Table  3  reveals  that  the  number  of  required  iterations  depends  mainly  on  the 
side-to-wavelength  ratio  of  the  cube,  and  hardly  at  all  on  the  number  of  patches  per 
square  wavelength  or,  equivalently,  hardly  at  all  on  the  number  of  unknowns  N  for  a 
fixed  sA  (assuming  a  reasonable  minimum  number  of  patches  per  square  wavelength  are 
used).  This  independence  of  the  number  of  conjugate  gradient  iterations  on  the  cell 
density  has  also  been  observed  in  the  solution  to  two-dimensional  scattering 
problems.3- 1 41  Since  the  CPU  time  is  proportional  to  the  total  number  of  patches,  the  CPU 
time  for  the  conjugate  gradient  solution  is  minimized  by  choosing  the  least  number  of 
patches  per  square  wavelength  sufficient  for  the  desired  solution  accuracy  (approxi¬ 
mately  25  patches  per  square  wavelength  in  our  case). 

The  number  of  iterations  Ii  required  for  convergence  with  the  conjugate  gradient 
method  as  a  function  of  the  number  of  unknowns  N,  when  using  the  minimum  patch  density 


Wheeler  III,  J.E.  and  Wilton,  D.R.  (1988)  Comparison  of  convergence  rates  of  the 
conjugate  gradient  methoJ  applied  to  various  integral  equation  formulations,  Diges!  of 
the  IEEE  AP-S  Symposium,  Syracuse,  NY,  pp.  229-232. 

■» 

Peterson,  A.F.  and  Mittra,  R.  (1986)  Convergence  of  the  conjugate  gradient  method 
when  applied  to  matrix  equations  representing  electromagnetic  scattering  problems,  IEEE 
Trans  Antennas  Propagat.  AP-34:  1447-1454. 

"  Peterson,  A.F.,  Smith,  C.F.,  and  Mittra,  R.  (1988)  Eigenvalues  of  the  moment-method 
matrix  and  their  effect  on  the  convergence  of  the  conjugate  gradient  alrorithm,  IEEE 
Trans  Antennas  Propagat  36:1177-1179. 


of  about  25  per  square  wavelength  [N  =75(sA)2],  appears  from  Table  3  to  approach  a 
logarithmic  function  of  N  as  N  gets  large;  specifically 


I  =  33  ln(.02N)  =  66  In(  1 .25  sA).  (33) 

Substituting  I  for  I  in  Eq.  (29)  [divided  by  .75  since  Eq.  (29)  gives  a  predicted 
value  that  is  about  75  percent  of  the  actual  CPU  run  time]  gives 


T  *  8.5  x  10  6  N2ln(.02N)  =  .l(sA)4  in(sA) 
co 


(34) 


minutes,  as  an  estimate  of  the  CPU  time  required  to  solve  for  scattering  from  large 
cubes  by  the  conjugate  gradient  method  on  our  32-bit,  1.3  megaflop  Linpack- 
performance-rated  computer.  (Interestingly,  Catedra  et  al.34  also  found  a  CPU 
time  dependence  proportional  to  the  right  side  of  Eq.  (34)  when  solving  3-D  scattering 
problems  using  the  conjugate  gradient  fast  Fourier  transform  method  applied  to  a  volume 
electric- field  integral  equation.)  Because  the  logarithmic  function  is  so  slowly  vary¬ 
ing,  Eq.  (34)  implies  that  the  CPU  time  for  solving  full-bodied,  3-D,  multi-wavelength 
scatterers  with  well-behaved  surface  integral  equations  increases  roughly  as  the  fourth 
power  of  the  electrical  size  of  the  scatterer. 

In  Figure  4  the  conjugate  gradient  and  Gaussian  elimination  CPU  times  vs  sA  for 
scattering  from  the  cube  are  plotted  from  Eqs.  (34)  and  (30)  [with  N  =  75(sA)2]  by  the 
solid  and  dashed  lines,  respectively.  Even  though  two-fold  symmetry  of  the  cube  has 
been  utilized  to  reduce  the  number  of  unknowns  N  by  the  factor  of  4,  Figure  4  confirms 
that  Gaussian  elimination  CPU  time  becomes  prohibitive  for  cubes  larger  than  a  few 
wavelengths  across,  and  that  conjugate  gradient  iteration  allows  one  to  determine 
scattering  from  considerably  larger  bodies. 


34  Catedra,  M.F.,  Gago,  E.,  and  Nuno,  L.  (1989)  A  numerical  scheme  to  obtain  the  RCS 
of  three-dimensional  bodies  of  resonant  size  using  the  conjugate  gradient  method  and 
the  fast  Fourier  transform,  IEEE  Trans.  Antennas  and  Propagat.,  37:  528-537. 
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Figure  4.  CPU  Time  Versus  Side-Length  to  Wavelength 
Ratio  for  Dual-Surface  Magnetic-Field  Integral  Equation 
Solution  for  Scattering  from  a  Perfectly  Conducting 
Cube.  The  number  of  unknowns  N  is  given  by  75  (sA)  , 
since  two-fold  symmetry  of  the  cube  was  used  to  reduce 
the  number  of  unknowns  by  a  factor  of  four  and  the  fixed 
patch  density  is  25  per  square  wavelength.  The  conjugate 
gradient  and  Gaussian  elimination  times  shown  here  would 
be  reducetl  by  factors  of  about  two  and  five,  respectively, 
if  the  coefficient  matrices  could  be  stored  in  central 
memory  rather  than  in  direct  access  files,  that  is, 
if  the  CPU  time  were  computed  from  Eqs.  (32)  and  (31) 
instead  of  Eqs.  (34)  and  (30). 
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We  emphasize  that  the  formula  (34)  for  the  conjugate  gradient  CPU  time  as  a 
function  of  the  number  of  unknowns  and  electrical  size  of  the  scatterer  is  an 
approximation  obtained  by  solving  for  scattering  from  the  perfectly  conducting  cube 
using  the  magnetic-field  dual-surface  integral  equation.  The  fonnula  holds  for  a  patch 
density  of  about  25  patches  pet  square  wavelength  and  a  normalized  residual  of  10 6. 

Since  the  number  of  iterations  is  nearly  independent  of  patch  density,  higher  patch 
densities  will  increase  the  CPU  time  proportionately.  The  CPU  time  will  decrease  if  the 
normalized  residual  is  chosen  greater  than  10  6;  in  particular,  we  found  that  the  number 
of  iterations  and  thus  CPU  time  halved  when  the  normalized  residual  was  increased  from 
106  to  a  value  between  10 3  and  10  2.  As  mentioned  in  Section  2,  the  number  of 
iterations  and  CPU  time  wdl  also  vary  somewhat  with  the  chosen  values  of  the  parameters 
a  and  8  in  the  dual-surface  integral  equation,  but  this  variation  was  not  large  for 
values  of  a  between  ±1.5  i  and  8  between  X/8  and  3 A/8. 

We  also  applied  the  magnetic-field  dual-surface  integral  equation  to  rectangular 
boxes  with  side-length  ratios  that  differed  considerably  from  the  value  of  unity  for  the 
cube.  For  some  rectangular  boxes,  the  required  number  of  iterations  and  CPU  time  were 
appreciably  larger  than  the  values  predicted  by  Eqs.  (33)  and  (34)  for  a  cube  of  the  same 
surface  area,  but  T  in  Eq.  (34)  was  never  larger  than  8.5  x  10 6  N5/2.  Although  the 
incident  plane  wave  always  propagated  normally  to  the  xy  face  of  the  rectangular  boxes 
(broadside  incidence),  it  is  unlikely  that  the  N-dependence  of  T^  in  Eq.  (34)  would 
change  dramatically  with  the  direction  of  the  incident  plane  wave,  because  the 
formulation  took  advantage  of  the  xy  symmetry  that  results  from  the  broadside  incidence 
to  reduce  the  number  of  unknowns  N  in  the  coefficient  matrix  by  a  factor  of  4. 

Finally,  in  hopes  of  reducing  computer  time  further,  we  experimented  with  three 
variations  of  the  conventional  conjugate  gradient  method,  namely  the  "biconjugate" 
gradient  method,  the  "augmented"  conjugate  gradient  method,  and  the  "modified"  conjugate 
gradient  method.  We  found  that  for  the  three  dimensional,  multi-wavelength  problem 

solved  with  surface  integral  equations,  these  tliree  variations  converged  more  slowly 
than  the  conventional  conjugate  gradient  method,  regardless  of  the  initial  guess,  or 
whether  they  were  used  alone  or  in  conjunction  with  the  conventional  conjugate  gradient 
method.3 


35  Sarkar,  T.K.  (1987)  On  the  application  of  the  generalized  biconjugate  gradient 
method,  J  Electromagnetic  Waves  and  Applications,  1(3):  223-242. 
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